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1. Introduction 

In these lectures, I will briefly review some of the Quantum Monte Carlo (QMC) 
methods that have been used to calculate properties of the "electron gas" and review 
properties that have been computed with these methods. (Note that this article is a 
mosaic of previously published works on this subject with very little genuinely new 
results. It is hoped that this collection will prove useful.) For the electron gas, QMC has 
been very useful in providing exact results, or at least exact constraints, on its properties, 
though with some remaining approximations. Fermion statistics remain a challenge to 
the practitioner of quantum simulation techniques. Nonetheless, the results are in many 
cases more accurate than those from the other quantum many-body methods, particularly 
for the case of strong correlation. There has been significant progress in the last decade 
in improving the methods. We will describe these improved methods. 

The ground state properties of the electron gas are entirely determined by the density 
parameter r s — a/a^ where in 3D, 4irpa 3 /3 — 1 and ao is the bohr radius, possibly 
changed from its vacuum value by band effects. In effective Rydbergs, the Hamiltonian 
is : 



Note that the kinetic energy scales as l/r 2 s and the potential energy scales as l/r s so 
that for small r s (high electronic density), the kinetic energy dominates, and the electrons 
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behave like ideal gas; in the limit of large r s , the potential energy dominates and the 
electrons crystallize into a Wigner crystal[l]. 

First few words on notation. I will always assume that the system is a non-relativistic 
collection of TV particles with A = h 2 /2m. I will stick to the first quantized notation 
in the canonical ensemble. A boson wave function is then totally symmetrical under 
particle exchange and a fermion function is antisymmetrical. The symbol R refers to the 
37V vector of particle coordinates, a to the N spin coordinates, and (r^, a) to the 3 spatial 
and 1 spin coordinate of particle i. The permutation operator acting on particle labels 
is denoted PR. Sometimes I will refer to the exact eigenfunctions and eigenvalues of the 
Hamiltonian: {(j) a {R),E a ). A known (computable) trial wave function will be denoted 
^{R). The symbol J will imply an integral over the full 37V dimensional configuration 
space of the particles, and possibly a sum over spin variables. 

We first discuss the QMC methods at T=0, then the obtained results, then the finite 
temperature methods and corresponding results. We restrict ourselves to the homoge- 
neous electron gas. There are many applications to inhomogeneous systems [3] but we do 
not have space to consider them. 

2. Variational Monte Carlo 

2T. Random Walks. - Monte Carlo methods for many-body systems are exclusively 
examples of Markov processes or random walks. Let us briefly review their basic concepts. 
Let S be a state space and let (so, si, . . .) be a random walk in that space. The choice 
of the state space will vary depending on the method, but for now, let s t represent the 
37V coordinates of all the particles. In the simplest Markov process, there is a constant 
transition probability for generating state b given that the walk is presently in state a, 
which we will denote Pf, a . The transition probability, or moving rule, generates the walk. 
Under very general conditions [2], the asymptotic probability distribution of the walk 
converges exponentially fast to a unique distribution: 



(2) V(s n ) ->II(a n ). 

In projector Monte Carlo, the transition rules are set up so that the asymptotic 
population is the ground state wave function for a given Hamiltonian. In the Metropolis 
[4, 5] rejection method, moving the particles is a two step process; first one samples a 
trial position, (state a'), from a transition probability T a > a , then this trial move is either 
accepted (i.e. b = a') or rejected (i.e. b — a) with a probability given by: 

(3) min 

The acceptance probability has been chosen to satisfy the detailed balance relation which 
implies that its asymptotic probability will converge to II a independent of the transition 
rules T aa i . The rate of convergence or the efficiency of the walk in sampling state space 
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is determined by the transition rules and the distribution to be sampled. The rejection 
method is appropriate when one wants to sample a known, computable function. If one 
had an exact analytic expression for the many-body wave function, it would then be 
straightforward to use this method to determine quantum expectation values in that 
state. However, such is not the case, and one is forced to resort to either more compli- 
cated, or more approximate methods. 

2'2. The Variational Method. - The variational Monte Carlo method was first used by 
McMillan [6] to calculate the ground state properties of liquid 4 He and then generalized 
to fermion systems by Ceperley et. al[7]. The variational theorem says that for ^ a 
proper trial function, the variational energy of the trial function is an upper bound to 
the exact ground state energy: 

u) F _ aH) >F 

(4) Ev s**(R)m) 

The trial function must satisfy the following conditions: 

1. has the proper symmetry: ^(i?) = (— l) Pl $>(PR) for fermions and the right 
behavior at the periodic boundaries. 

2. is well defined everywhere which means that both \P and V'J must be continuous 
wherever the potential is finite. 

3. The integrals J* 2 , j^ 2 H^, and j(^H) 2 should exist. The last integral is only 

required to exist for a Monte Carlo evaluation of the integrals. If it does not exist 
the statistical error of the energy will be infinite. 

In the continuum, it is important to show analytically that properties 2 and 3 hold 
everywhere, particularly at the edges of the periodic box and when two particles approach 
each other. Otherwise either the upper bound property is not guaranteed or the Monte 
Carlo error estimates are not valid. 

The variational method is then quite simple. Use the Metropolis algorithm to sample 
the square of the modulus of the wave function: 

(5) n<«)= l*™> 2 



/ 1 *(fi) 



(Note that this may also include spin degrees of freedom.) Then the variational energy 
is simply the average value of the local energy over this distribution, 

(6) E v = J U{R)E L {R) = (E L (R)) 

where the local energy of is defined as: 

(7) E L (R) = 
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Variational Monte Carlo(VMC) has a very important zero variance property : as the trial 
function approaches an exact eigenfunction, the local energy approaches the 

eigenvalue everywhere. El(R) — ► E a . and the Monte Carlo estimate of the variational 
energy converges more rapidly with the number of steps in the random walk. Of course, 
in this limit the upper bound is also converging to the true energy. It is because of the 
zero variance property that Quantum Monte Carlo calculations of energies can be much 
more precise than Monte Carlo calculations of classical systems. Fluctuations (statistical 
errors) are only due to inaccuracies in the trial function. 

2'3. The Pair Product Trial Function. - The pair product trial function is the simplest 
generalization of the Slater determinant and the ubiquitous form for the trial function 
in variational Monte Carlo: 

(8) *(R,<r) = ea;p[-5^u(r tf )]det[0fc(»"i,<ri], 

i<j 

where 9k(r,<r) is the k th spin-orbital and ^(r,, <7j), the Slater matrix, and u(r) is the 
"pseudopotential" or pair correlation factor. This function also goes by the name of a 
Jastrow[8] wave function, although Bijl[9] much earlier described the motivation for its 
use in liquid 4 He. Closely related forms are the Gutzwiller function for a lattice, or the 
Laughlin function in the fractional quantum hall effect. Both u{r) and 9k(r,a) arc, in 
principle, determined by minimizing the variational energy. 

2'4. Details. - I will only mention a few details concerning VMC. First, how do the 
particles move? In the continuum it best to move the electrons one at a time, by adding 
a random vector to the electron's coordinate, where the vector is cither uniform inside 
of a cube centered about the old coordinate, or is a normally distributed random vector, 
as we will discuss below with DMC. Assuming the first kind of move for the i th particle, 
the trial move is accepted with probability: 

(9) | I 2 - exp^^Wrl Ti) u{n - r,))] ^WWki | 2 , 

jjti k 

where the matrix, C, is the transposed inverse to the Slater matrix. The evaluation of a 
general determinant takes 0(N 3 ) operations. The evaluation of the fcrmion part of the 
acceptance ratio will take O(N) operations if C is kept up to date. If a move is accepted, 
C needs to be updated [7] which takes 0(N 2 ) operations. Hence, to attempt moves for 
all N particles (a pass) takes 0(N 3 ) operations. 

The local energy, needed to evaluate the variational energy, is calculated by applying 
the Hamiltonian to the trial function: 



(10) 



E L (R) = V(R) + A]T[V 2 L/ -^%(n)C ki - G 2 ] 

i k 
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where Gi = — ViU + ^2 k Vi6k(ri)Cki and U = X) w ( r ij)- Thus, the inverse matrix is 
also needed to determine the local energy. Very often the orbitals are taken to be exact 
solutions to some model problem, in which case the term, Vf^ri), will simplify. 

2'5. Optimization of Trial Functions. - Optimization of the parameters in a trial func- 
tion is crucial for the success of the variational method and important for the Projector 
Monte Carlo method. There are several possibilities of the quantity to optimize. 

• The variational energy: Ey. Clearly one minimizes Ey if the object of the cal- 
culation is to find the least upper bound to the energy. There are some general 
arguments suggesting that the trial function with the lowest variational energy will 
maximize the efficiency of Projector Monte Carlo [10]. 

• The dispersion of the local energy of the variance : J[(H — Ey) 1 ^] 2 ■ If we assume 
that every step on a QMC calculation is statistically uncorrelated with the others, 
the dispersion is proportional to the variance of the calculation. The minimization 
of the dispersion is statistically more robust than the variational energy because it 
is a positive definite quantity. 

• The overlap with the exact wave function: / ^f(j). This is equivalent to finding 
the trial function which is closest to the exact wave function in the least squares 
sense. This is the preferred quantity to optimize if you want to calculate correlation 
functions, not just ground state energies. Optimization of the overlap will involve 
a Projector Monte Carlo calculation to determine it, a more computer intensive 
step. 

We will now review the analytic properties of the optimal pseudopotential. Let us 
assume that the spin-orbits come from an exact solution for some model potential. Con- 
sider bringing 2 particles together and let us examine the dominant terms in the local 
energy. In a good trial function, the singularities in the kinetic energy must cancel those 
in the potential. The singular parts of local energy will have the form: 

(11) E L (R) = v(r) + 2AV 2 u(r) - 2A(Vw(r)) 2 + . . . 

where r is the distance separating the particles. An intuitive result emerges: e~ u ^ will 
equal the solution to the 2-body Schrodinger equation. For the Coulomb potential this 
equation gives the "cusp condition" , the value of the derivative of the pseudopotential 
at the origin: 



e 2 
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Even more important is the large r behavior of the optimal u(r) where a description in 
terms of collective coordinates (the plasmons) is appropriate. This is important because 
the long-wavelength modes arc important for the low energy response properties and 
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they are also the slowest to converge in QMC. In k-space, the variational energy can be 
written as: 

(13) E v = E F + J2( v k-\k 2 u k )(S k -l) 

k 

where Ep is the fermion energy in the absence of correlation, Vk and Uk are the fourier 
transforms of v(r) and u(r) and S k is the static structure factor for a given u(r). Min- 
imizing Ey with respect to Uk and making the RPA assumption of how Sk depends on 
Uk : S^ = SQ k + 2pUk where p is the particle density and Sok is the structure factor for 
uncorrelated fermions, we obtain [11] the "Gaskell" or RPA form[12] at long wavelengths: 



(14) 2pu k 



1 + jj_ + 2 ^fcji/ 2 
Sok Sctk Afc 2 



Here the fourier transform of the pair potential is v k — 2ir(D — ^e 2 /? 13-1 . This gives the 
following behavior of the 2DEG and 3DEG pseudopotential: lim r ^oo u(r) cx r^- 0-1 )/ 2 . 
Not only the power, but also the coefficient is exactly correct. This Gaskell form is 
optimal at both long distances and at short distances, but not necessarily in between. 
For the 2DEG and 3DEG, this parameter-free form can hardly be improved by assuming 
an arbitrary function and is basically correct even in the strongly correlated regime, 
though the RPA assumption on which it is based, is not. 

This raises a very important point which we will not have space to go into. The 
Coulomb potential is long-ranged, so that treatment of a finite system in periodic bound- 
ary conditions requires an assumption about how to handle the part of the potential 
sticking out of the box. For high accuracy results and physically correct properties in 
the long wavelength limit, the Ewald image method[13, 11] is needed to represent the 
correct long-range behavior of the potential and the trial function. The optimal pseu- 
dopotcntials are always long-ranged in the sense that correlation will extend beyond the 
simulation box. The ground state energy is not very much affected by this tail in the 
wave function, but response function, such as the dielectric function or static structure 
factor are crucially dependent on using the correct long-range properties. In order to 
maintain the upper bound property, the correlation function must be properly periodic 
in the simulation cell. See rcf. [14] for an efficient way to treat numerically any long-range 
function within periodic boundary conditions. 

At intermediate distances or for highly correlated or complex systems, a purely Monte 
Carlo optimization method may be needed. The simplest such method consists of running 
independent VMC runs with a variety of different variational parameters, fitting the 
results energies with a quadratic form, doing more calculations at the predicted minimum, 
etc. , until convergence in parameter space is attained. The difficulty is that close to 
the minimum the independent statistical errors will mask the variation with respect to 
the trial function parameters. The derivative of the variational energy with respect to 
trial function parameters is very poorly calculated. It is difficult to optimize by hand 
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functions involving more than 3 variational parameters. A correlated sampling method, 
known as reweighting[15, 7] solves this problem. 

2'6. The spin-orbits: Crystal vs. Liquid. - The description of the trial function in the 
previous section did not discuss the spin-orbits that enter into the slater determinant nor 
how we will distinguish between the liquid and crystal. In a homogenous liquid, plane 
waves are a self-consistent solution to the Hartree-Fock equations with the occupied 
states chosen to have an energy less than the Fermi level. It is known, however, that this 
solution is unstable to broken symmetry states, in particular to ferromagnetic solutions 
and to spin density waves. Attempts to use such orbitals in conjunction with a Jastrow 
factor have generally led to higher energies. We will discuss the ferromagnetic case where 
one does find a deviation from the non-interacting spin-orbits. 

At very strong correlations, Wigner showed the system must crystallize. However, 
the plane wave determinant even with a Jastrow correlation factor, does not provide a 
good description of the crystal. At the VMC level, the broken spatial symmetry of the 
crystal needs to be put in by hand. That is typically done by making spin-orbitals which 
are localized on lattice sites. It has been found the spherical Gaussian work as well as 
more complicated functions for the 2DEG and 3DEG: 



where C is a variational parameter and Zj is a set of lattice sites. To minimize the 
potential energy these are usually taken as a triangular lattice in 2D and a bcc lattice in 
3D. Since the bcc lattice is bipartite, one can then place the two species of spin on the 
two sublattices. The triangular lattice cannot be handled this way. We will discuss the 
spin properties of the Wigner crystal at the end of these notes. 

It is found that these spin-orbits are very successful [32] for the perfect crystal, even 
though the broken symmetry is put in by hand. Note that the long range properties 
of the optimal Jastrow factor are slightly modified[ll]. "Shadow" trial functions, where 
one introduces a set of shadow particles [16], can spontaneously break the symmetry, and 
provide a better description of defective crystals, or near melting, however, treatment of 
Fermi statistics with those functions is numerically difficult. 

2 7. Beyond the pair-product trial function. - There is important progress in finding 
trial functions significantly more accurate than the pair product, Slater-Jastrow function. 

• One can replace the spin-orbitals with pairing functions. For example, if the par- 
ticles are paired in a spin singlet state one obtains a determinant of the form: 



(15) 



0i(r) = exp(-C(r - Z,) 2 ) 



(16) 



detain ],r } 



I)] 



or if the particles are paired in a spin triplet: 



N/2 



(17) 
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where A is an antisymmetrizer. Bouchaud and Lhuiller [17] have pointed out the 
this object is a pfaffian and hence one can use the theorem that the square of a 
pfaffian is a determinant to sample it with VMC. To date, no one has succeeded in 
getting lower energy with this trial function than the SJ energy. 

• The dominant term missing in the trial function for bosonic system is a three body 
term with the functional form of a squared force: 

(18) W) = -£E£( r «)r«] 2 - 

» 3 

The form, a square, makes it particularly rapid to compute. It is not much slower 
than a 2 body function. This term is particularly important at strong correlation, 
at large r s [18, 19]. 

• The dominant correction for the pair-product trial function is to modify the spin- 
orbitals to include backflow correlations. The particle coordinates in the Slater 
determinants become 'quasi-particle' coordinates: 

(19) det[6 k (xi,<Ti], 

where the "quasi-particle" coordinates are defined by: x, = + ^2ijf){ r ij)^io- 
Backflow is needed to satisfy local current conservation. The r](r) function is similar 
to the "pscudopotcntial" however it decays more quickly, as r -3 in 3D. A good 
approximation to its form is discussed in [21]. Because this affects the orbitals, 
it also affects the nodes (or phases) of the many-body trial function. Thus it is 
extremely important in conjunction with the fixed-node approach. 

Table 1 gives VMC energies for the 3DEG with SJ trial function and with improved 
3 body backflow. There has been recent progress[21] in finding analytic properties of the 
backflow functions. Such analytic forms are quite accurate for r s < 20 and are shown in 
figs. 1-2. 

2'8. Finite size effects. - The dependence of the energy on the number of electrons 
is the largest systematic error in quantum Monte Carlo calculations but fortunately it 
can be minimized by several methods. The largest effect in fermion systems comes from 
the kinetic energy, arising from the erratic filling of the shells in k-spacc. If one uses 
periodic boundary conditions for the wavefunction the resulting energy per electron for 
free particles is shown in fig. 3 

Though PBC with Fermi Liquid Theory (FLT) corrections[ll] are adequate for unpo- 
larized and fully polarized systems, the precision is limited for intermediate polarizations. 
To estimate finite size effects within FLT, one must perform accurate DMC simulations 
for widely varying system sizes. For example, in Ortiz et al.[22] calculation of the elec- 
tron gas, the simulation size varied from 725 < N < 1450. Within DMC it is very 
time-consuming to ensure uniform accuracy independent of particle number, so that 
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Fig. 1. - The change in the quasiparticle coordinate rr/(r) (analytic backflow) caused by an 
electron a distance r away in the 3D electron gas. Graphed is only the short range part of r\ 
with iV = 54. The four figures are for r s = 1, 5, 10, 20 from the bottom to the top of the figure. 
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wavefunction 


E v 


a 


Edmc 


1 


SJ 


1.0669 (6) 


1.15 (2) 


1.0619 (4) 




BF-A 


1.0611 (2) 


0.029 (1) 


1.0597 (1) 


5 


SJ 


-0.15558 (7) 


0.0023(1) 


-0.15734 (3) 




BF-A 


-0.15762 (1) 


0.00061 (1) 


-0.15810 (1) 


10 


SJ 


-0.10745 (2) 


0.00039 (.5) 


-0.10849 (2) 




BF-A 


-0.10843 (2) 


0.00017 (1) 


-0.10888 (1) 


20 


SJ 


-0.06333 (1) 


0.000064 (1) 


-0.06388 (1) 




BF-A 


-0.06372 (2) 


0.000045 (2) 


-0.06408 (1) 



Table I. - Energies and variances for the 3D electron gas with N = 54 unpolarized electrons 
in Rydbergs/ electron. SJ means a Slater determinant of plane waves times an optimized Jas- 
trow factor. BF-A are the results using the RPA Jastrow, together with the analytical backflow 
formula, o 2 is the variance of the local energy per electron. 
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Fig. 2. - The three-body contribution to the logarithm of the wavefunction due to three electrons 
in the 3D electron gas. This is just the short range parts of w(r) for N = 54. The solid line, 
for r s = 1, is close to zero (maximum magnitude of 3 x 1CP 4 ). The dotted line and dashed lines 
are for r 3 — 5, 10. 



one typically determines size effects within VMC, using the more approximate SJ trial 
functions. 

Within periodic boundary conditions (PBC), the phase picked up by the wavefunction 
as a particle makes a circuit across the unit cell, is arbitrary. General boundary conditions 
are: 

(20) *(n + L,r 2 ,...) = e i9 *(ri,r 2) ...) 

where L is a lattice vector of the superccll. If the twist angle 9 is averaged over, most 
single-particle finite-size effects arising from shell effects in filling the plane wave orbitals, 
are eliminated [23] as can be seen in the figure. If one uses grand canonical simulations, 
all single particle kinetic finite size effects arc eliminated. TA boundary conditions allow 
a much better way to estimate energies in the thermodynamic limit of partially polarized 
fcrmi liquids since the number of electrons can be held fixed as the spin polarization 
varies. 
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Fig. 3. - Relative error of the energy versus number of particles with PBC and TABC in 2D 
and 3D. The points shown are only those where the relative error has a local maximum. Curves 
are shown only for N < 100. 



The potential effects are more difficult since they are due to pair interaction effects. 
To correct this, one fits the energies versus N using the expansion: 

(2D E N = E 00+ a ± + ^ + .... 

Results for several different sizes are necessary to use this formula, though we have made 
some recent progress with analytic corrections within the GCE. 

2'9. Problems with Variational Methods. - The variational method is very powerful, 
and intuitively pleasing. One posits a form of the trial function and then obtains an upper 
bound. In contrast to other theoretical methods, no further essential approximations 
need to be made and there are no restrictions on the trial function except that it be 
computable in a reasonable amount of time. To be sure, the numerical work has to be 
done very carefully which means that convergence of the random walk has to be tested 
and dependence on system size needs to be understood. To motivate the methods to 
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be described in the next section, let me list some of the intrinsic problems with the 
variational method. 

• The variational method favors simple states over more complicated states. One 
of the main uses of simulations is to determine when and if a zero-temperature 
phase transition will occur. As an example, consider the liquid-solid transition at 
zero temperature. The solid wave function is simpler than the liquid wave function 
because in the solid the particles are localized so that the phase space that the 
atoms explore is much reduced. This means that if you compare liquid and solid 
variational energies for the same type of trial function, (e.g. a pair product form) 
the solid energy will be closer to the exact result than the liquid and hence the 
transition density will be systematically lower than the experimental value. VMC 
calculations in the e-gas show this effect [11]. Another example; the wave function 
for fully polarized electron gas is simpler than for unpolarized system so that the 
spin susceptibility computed at the pair product level has the wrong sign. 

• The optimization of trial functions for many-body systems is very time consuming, 
particularly for complex trial functions. This allows an element of human bias; the 
optimization is stopped when the expected result is obtained. 

• The variational energy is insensitive to long range order. The energy is dominated 
by the local order (nearest neighbor correlation functions). If one is trying to 
compare the variational energy of a trial functions with and without long range 
order, it is extremely important that both functions have the same short-range 
flexibility and both trial functions are equally optimized locally. Only if this is 
done, can one have any hope of saying anything about the long range order. The 
error in the variational energy is second order in the trial function, while any other 
property will be first order. Thus variational energies can be quite accurate while 
correlation function completely incorrect. 

• You almost always get out what is put in. Suppose the spin-orbitals have a Fermi 
surface. Then the momentum distribution of the pair product trial function will 
also have a Fermi surface although it will be renormalized. This does not imply that 
the true wave function has a sharp Fermi surface. Only if localized spin-orbitals 
are used will a gap appear. 

3. Diffusion Monte Carlo 

Now I will turn to a more powerful method than VMC where a function of the 
Hamiltonian projects out the the ground state, hence the name, projector Monte Carlo. 
For simplicity I will only discuss Diffusion Monte Carlo although most of what I say 
carries over immediately to other projectors. A sequence of trial functions is defined by 
applying the projector, G(R,R'): 



(22) 
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with the initial condition: ipo(R) = ^(R)- The effect on the trial function of the Green's 
function is seen by expanding the trial function in the exact eigcnfunctions of the Hamil- 
tonian: 

(23) MR) = 52 MR) < > e - n ^ E "- E -l 

a 

The Green's function will project out the state of lowest energy having a non-zero overlap 
with the initial trial function. 

(24) lim MR) = MR) < > e- nT{E °- ET \ 

n^oo 

The role of the trial energy is to keep the overall normalization of ip n fixed, which implies 
Et ~ E - The time step, r, controls the rate of convergence to the ground state. It is 
controlled by the needed to determine the projecting function. 

For many-body systems, the trial function and the Green's function are sampled. For 
the moment let us discuss the case where 'J is non-negative, the boson case. In the limit 
that the time step approaches zero, a coordinate space representation of the Green's 
function is: 

(25) <R\e-< n - E ^R> ^{^Xry^l^-^^e-^^-^ + 0(r 2 ) 

The iteration equation, eq. 25, has a simple interpretation in terms of branching ran- 
dom walks since the first factor is the Green's function for diffusion and the second is 
multiplication of the distribution by a positive scalar. 

An ensemble of configurations is constructed with a Metropolis sampling procedure for 
^f(R)- This is called the zeroth generation, i.e. n = 0, and the number of configurations is 
the population of the zeroth generation, P . Points in the next generation are constructed 
by sampling the Gaussian distribution in cq. (25) and then branching. The number of 
copies of R' in the next generation is the integer part of u + exp[— t(V(R) — Et)] where 
u is a uniform random number in (0, 1). If the potential energy is less than the ground 
state energy, duplicate copies of the configuration are generated. In future generations, 
these walks propagate independently of each other. In places of high potential energy, 
random walks are terminated. 

The above procedure, depicted in fig. (2), is a Markov process where the state of 
the walk in the nth generation is {P„; R\, R2, ■ ■ ■ Rp n }- Hence it has a unique stationary 
distribution, constructed to be the ground state wave function. The number of walk- 
ers fluctuates from step to step. The trial energy, Et, must be adjusted to keep the 
population within computationally acceptable limits. 

31. Importance Sampling. - The above scheme, first suggested by Fermi, was actu- 
ally tried out in the first days of computing some fifty years ago [24]. But it fails on 
many-body systems because the potential is unbounded. For example, a coulomb po- 
tential can go to both positive and negative infinity. Even with a bounded potential 



14 



D. M. Ceperley 



the method becomes very inefficient as the number of particles increases. But there is a 
very simple cure discovered by Kalos[25] for GFMC but equally applicable to any pro- 
jector method. Importance sampling multiplies the underlying probability distribution 
by a known, approximate solution. Multiply eq. (22 ) by "f, the trial function, and let 
f n (R) = V(R)MR)- Then: 

(26) f n+1 = *e-^-^)Vn = / dR'G(R, R')f n (R') 

where G(R,R') = ^/~ 1 e~ T ^ H ~ ET ^ is the importance-sampled Green's function and the 
initial conditions arc fo{R) = *f? 2 (R). It is easily shown by differentiating G with respect 
to r that it satisfies the evolution equation: 

(27) _ 9G(R^Ro;r) = _ ^ A . v . [v .£ + 2 £ v . In(*(i2))] + [E L (R) - E T ]G 

i 

where El is the local-energy defined in cq. (7). The first three terms on the right-hand 
side correspond to diffusion, drifting and branching. As the trial function approaches 
the exact eigenfunction, the branching factor approaches unity; thus a sufficiently good 
trial function can control the branching. 

The importance sampled DMC algorithm[26, 51] is: 

1. The ensemble is initialized with a VMC sample from \^>(R)\ 2 (R) 

2. The points in the configuration are advanced in time as: 

(28) R n+1 = Rn + X + ArV ln(\<f (R n )\ 2 ) . 

where x is a normally distributed random vector with variance 2At and zero mean. 

3. The number of copies of each configuration is the integer part of 

(29) exp(-T(E L (R n ) - E T )) + u 
where u is a uniformly distributed random number in (0, 1). 

4. The energy is calculated as the average value of the local energy: E n = (_E^(i?„)). 

5. The trial energy is periodically adjusted to keep the population stable. 

6. To obtain ground state expectations of quantities other than the energy, (e.g. the 
potential energy) one must correct the average over the DMC walk, the so-called 
"mixed estimator", V m i X —< ^ol^l* > with using the variational estimator [15]. 

(30) < (f> \V\(f> >= 2 < (f>o\V\^ > - < > +O([0 O - *] 2 ) 

If the mixed estimator equals the variational estimator then the trial function has 
maximum overlap with the ground state. 
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3'2. The Fixed-Node Method. - We have not discussed at all the problem posed by 
fermi statistics for projector Monte Carlo. Consider the difficulty in implementing the 
non-importance sampled algorithm: the initial condition is not a probability distribution 
since a fermi on trial function will have positive and negative pieces. Hence, we must use 
the initial sign of the wave function as a weight for the random walk. That leads to an 
exact, but slowly converging algorithm, that we will discuss it in the next subsection. 

Now suppose we use an antisymmetric trial function in the importance-sampled al- 
gorithm. The initial distribution is positive, but the Green's function, G(R, R') can be 
negative if a step changes the sign of vp. Thereafter a minus sign will be attached to 
the walk which will lead to a growing statistical variance for all estimators. But there 
is a simple way to avoid the sign: forbid moves in which the sign of the trial function 
changes. This is the fixed-node (FN) approximation. 

In a diffusion process, forbidding node crossings puts a zero boundary condition on the 
evolution equation for the probability. This solves the wave equation with the boundary 
conditions that it vanish wherever the trial function vanishes. One can easily demonstrate 
that the resulting energy will be an upper bound to the exact ground state energy [27]; the 
best possible upper bound with the given boundary conditions. With the FN method, 
we do not necessarily have the exact fermion energy, but the results are much superior 
to those of VMC. No longer do we have to optimize two-body correlation factors, three- 
body terms etc., since the nodes of the trial function are unchanged by those terms. 
One is exactly solving the wave equation inside the fixed-nodal regions, but there is a 
mismatch of the derivative of the solution across the boundary. The nodes have an un- 
equal "fermion" pressure on the two sides unless the nodes are exact. Where comparison 
has been done between the VMC energy, the FN-DMC energy and the exact answer, 
one generally finds that the systematic error in the FN calculation is three to ten times 
smaller than it would be for a well optimized VMC energy. 

The nodes obviously play a very important role since, as we have seen, if the nodes 
were exactly known, the many-fermion system could be treated by Monte Carlo methods 
without approximation. The ground state wave function can be chosen real in the absence 
of magnetic fields; the nodes are the set of points where 4*{R) — 0. Since this is a single 
equation, the nodes are in general a 3N — 1 dimensional hypersurface. When any two 
particles with the same spin are at the same location the wave function vanishes. These 
coincident planes, with = rj are 3iV — 3 dimensional hypcrsurfaces. In 3D space they 
do not exhaust the nodes, but are a sort of scaffolding. The situation is very different in 
ID where the set of nodes is usually equal to the set of coincident hyperplanes. Fermions 
in ID arc equivalent to ID bosons with a no-exchange rule. 

Nodal volumes of ground state wave functions possess a tiling property [28]. To define 
this property first pick a point, Ro, which does not lie on the nodes. Consider the set 
of points which can be reached from i? by a continuous path with <j>(R) ^ 0. This is 
the volume in phase space accessible to a fixed-node random walk starting at Rq. Now 
consider mapping this volume with the permutation operator (only permute like spins), 
i.e. relabel the particles. The tiling theorem says that this procedure completely fills 
phase space, except, of course, for the nodes. Thus one does not have to worry about 
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where the random walk started; all starting places are equivalent. This theorem applies 
for any fcrmion wave function which is the ground state for some local Hamiltonian 
and it can be proved by a simple variational argument. Excited states, or arbitrary 
antisymmetric functions need not have the tiling property. More extensive discussion of 
fcrmion nodes and some pictures of cross-sections of free particle nodes are given in rcf . 
[28]. 

3'3. Fixed phase method. - We can generalize the fixed-node method to complex trial 
functions [43]. Let us take an arbitrary complex wave function and divide its logarithm 
into its real and imaginary parts: 4>{R) = exp(—M(R) + iS(R)). For a given fixed-phase, 
let us minimize the variational energy with respect to M. We are left with an effective 
Hamiltonian with an additional repulsive term in the potential energy coming from the 
phase of the trial function: 

N 

(31) AH = -\J2[VS(R)} 2 

i 

with S(R) the phase of the trial function. Then this equation can be solved with the 
bosonic DMC method described above. We do not solve the full Schrodinger equation 
since we have not optimized the energy with respect to the phase. That is we have an 
error in the continuity equation. In contrast to the method for real wave functions, there 
is no barrier to the diffusion (no nodes), just regions where the phase is rapidly varying, 
leading to a large (but not infinite) effective potential. 

3'4. Exact Fermion Methods. - As accurate as the FN method might be, it is still 
unsatisfactory since one does not know how the assumed nodal structure will affect the 
final result. One might guess that long-range properties, such as the existence or non- 
existence of a fcrmi surface will be determined by the assumed nodes. The FN algorithm 
is only improving the bosonic correlations of the trial function, not the fermion features. 
There are some fairly simple ways of improving on the FN method, but their use is 
limited. 

The transient estimate (TE) method calculates the ratio: 

(32) E TE (t) = J * rte 



/ ^ exp[-t(7i - E T )}^ 

Clearly the variational theorem applies: E TE (t) > E . Also the energy converges expo- 
nentially fast in t: 

(33) lim E TE {t) = E + 0(e- tE *) 

i— >oo 

where E g is the gap to the next excited state with the same quantum numbers as the 
fcrmion ground state. In a Fermi liquid this is the gap to the state with the same 
momentum, parity and spin and would be obtained by making 2 particle-hole excitations. 
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For a method to self-consistently find its own nodes, the walks must be able to go 
anywhere, and so the drift term in eq. (27) must not diverge at the nodes. Hence we 
must distinguish between the antisymmetric trial function that is used to calculate the 
energy, ^(i?), (this is always assumed to be our best variational function) and a strictly 
positive guide function, ^g(R), used to guide the walks. The guide function appears 
in the drift and branching terms of eq. (27) and will be assumed to be a reasonable 
boson ground state trial function while the trial function appears in eq. (32). The ^>g 
importance-sampled Green's function is: 

(34) G(R,R';t) = V G {R) < R\ e -^ n - E ^\R' > ^(R'), 

and we can rewrite eq. (32) as: 

/ a(R)E LT (R)G(R, R'\ t)a(R')^ G (R') 



(35) E TE {t) 



J a(R)G(R, R';t)a(R')^ 2 G (R') 



where a(R) = ^(R) /^ G (R) and E LT {R) is the local energy of In the limit, * G — ► |*|, 
a equals the sign of the trial function. 
The transient estimate algorithm is: 

1. Sample configuration R' from the square of the guide function with VMC. That 

corresponds to the rightmost factor in the above integrands. 

2. Record the initial sign of the walk, a(R'). 

3. Propagate the walk forward an amount of time, t with the Green's function, G(R, R'; t). 

If a branch occurs, each branch will count separately. 

4. The weight of the walk arriving at R is a(R)<r(R'). The energy at time t is computed 

as: 

,^ p m < [Elt(R) + E LT {R')]o{R)a{R>) > 
(36) EMt) = < a(R)a(Bf) > ' 

where the averages are over all random walks generated by this process. 

The weight of the walk is positive if the walk crosses an even number of nodes (or 
does not cross at all) and is negative if it crosses once or an odd number of times. The 
trial function nodes are displaced by an unequal diffusion of walks from one side or the 
other. 

The release node (RN) algorithm[27, 29] is an improvement on this TE method. 
Instead of projecting from the trial function, one begins the projection from the fixed- 
node solution. There are several advantages. First of all boson correlation within the 
fixed-nodes is already optimized, thus the projection time is only determined by the 
time to adjust the position of the nodes (of course this will indirectly affect the bosonic 
correlation). Second, one can directly calculate the difference between the exact result 
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and the fixed-node solution. It turns out that this is given by the local energy of walks 
as they cross the nodes. Thus the difference is obtained with more statistical accuracy 
than cither energy alone which allows the convergence to be carefully monitored. Finally 
the release node method can be conveniently integrated into a fixed-node program. The 
only modifications are to introduce a guide function, and to keep track of the energy as a 
function of time since nodal crossing. The Ceperley-Alser results[26] for the 3DEG were 
done with the RN method. For complex trial function one can perform released phase 
simulation [48], a straightforward generalization of release node. 

However there are serious problems with both the TE and RN method. Let us examine 
how the statistical error of the eq. (32) depends on the projection time. Note that 
the value of both the numerator and denominator arc asymptotically proportional to 
exp(— t(Ep — Et))- Thus to keep the normalization fixed our trial energy must be equal 
to Ep. But because the guide function allows the walks to cross the nodes, the population 
will increase as exp(— t(Es — Et)) where Eb is the boson energy. One can demonstrate 
that the signal-to-noise ratio vanishes exponentially fast. This is a general result. In any 
fermion scheme, as soon as negative weights are introduced the statistical error will grow 
as: 



(37) 



-t(E F — E B ) 



The behavior is physically easy to understand. Our estimator depends on finding differ- 
ences between random walks crossing an even or an odd number of times. As soon as 
there is substantial mixing, the difference becomes harder and harder to see. Note that 
the exponential growth rate depends on a total energy difference. This implies that the 
transient estimate algorithm is guaranteed to fail if N is sufficiently large; the statistical 
errors will be too large. Nonetheless reliable results have been obtained for systems of 
54 fcrmions. 

The convergence problem is actually a bit more subtle since the projection time, t, 
can be chosen to give approximately equal statistical errors and systematic errors coming 
from non-convergence of the projection. Taking these errors from eqs. (33,37) we find 
the total error will decrease as: 

where P is the total number of steps in the random walk. Only for bosons will rj = 1/2. 
Any excited state will converge at a slower rate. Note that r\ cx 1/N for a fermion system. 
Inverting this relation, we find that the computer time needed to achieve a given error 
will increase exponentially with N . 

There have been many attempts to "solve" the fermion sign problem. For example 
an obvious method is to try and pair positive and negative random walks in the TE 
method. This is difficult in many dimensions simply because the volume of phase space 
is so large than random walks rarely approach each other. There is some confusion about 
the nature of the "fermion" problem in the literature. Note that the TE and RN methods 
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do converge to the exact fermion energy. The fcrmion problem has to do with how long 
it takes to achieve a given error estimate, and, more precisely, how this scales with the 
number of fermions: the computational complexity of the calculation. 

3'5. Problems with Projection methods. - The projection method shares many of the 
same problems with the variational method. In fact, it is useful to think of the projection 
method as a "super- variational" method. In both the VMC and DMC methods there 
is a premium for good trial functions; that is the most straightforward way of making 
progress to solving the many-fcrmion problem. Our recent work[21] on finding more 
accurate analytic backflow-threebody shows that even for the well-studied homogeneous 
electron gas, significant improvements can still be made. Some general problems with 
projection methods: 

• The fixed-node result is guaranteed to be closer to the exact answer than the 
starting variational trial function. Since the FN algorithm automatically includes 
bosonic correlation, the results are much less likely to have the human bias than 
with VMC. There is also the possibility of new things coming out of the simulation. 
For example, one may observe a particular type of correlation completely absent 
from the trial function. Hence, it is always good to pay close attention to correlation 
functions computed by DMC since that it is a good way of learning what is missing 
in the trial function. But it is slower than VMC because the time step needs to be 
smaller. 

• Although the probability distribution does converge to the exact answer, in prac- 
tice, this does not always occur in any given calculation of a many-body systems. 
The situation is similar to that of a classical simulation near a phase boundary. 
Metastable states exist and can have a very long lifetime. However, with DMC 
the importance sampling always biases the result. If the trial function describes 
a localized solid, even after complete convergence, the correlation functions will 
show solid-like behavior. Careful observation will reveal liquid-like fluctuations in- 
dicating the presence of the other state. The ability to perform simulations in a 
metastable state is useful but the results must be interpreted with caution. 

• Importance sampling is only a partial cure to the unbounded fluctuations of the 
branching method. As N increases, sooner or later the branching becomes uncon- 
trollable. Most projector Monte Carlo calculations have fewer than several hundred 
fermions. Finite temperature Metropolis methods do not suffer from the problem 
of uncontrolled branching. Also, as system gets larger, it becomes more important 
to monitor convergence, since necessary projection times can become larger. One 
is doing a simulation in 4D space-time and as the spatial length L is scaled, the 
length in Euclidean time (i must be kept proportional. 

• Although the fixed-node approximation dramatically improves energies, other prop- 
erties, such as the momentum distribution, may not be improved. The projector 
methods can only calculate energies exactly. For all other properties one must 
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extrapolate out the effect of the importance sampling. This is a real problem if 
one is interested in obtaining asymptotic behavior of correlation functions. There 
are ways of getting around some of these problems such as forward walking, but 
none are totally satisfactory. The Path Integral finite temperature methods are 
superior to Projector Monte Carlo for calculating correlation functions. To explore 
the metal-insulator phase transition with FN-DMC, one must come up with a se- 
quence of nodes spanning the transition and use the upper bound property of the 
fixed-node approximation. Reptation Monte Carlo[31] allows another way around 
this problem of mixed-estimators, though the scaling versus problem size is still an 
issue. 

• Release node calculations only improve the nodes locally. If t is the release node 
projection time, then we can move the nodes a distance of at most y/6N\t. For 
this reason, it is important to build in exact long-range nodal properties. 

4. Ground state properties of the electron gas in 2 and 3 D 

Here I will briefly summarize some of the results coming from the use of variational 
and diffusion Quantum Monte Carlo methods for the electron gas. In the 3DEG the 
Wigner Crystal transition was estimated within VMC[11] as r s = 67. This was improved 
using DMC, both with fixed-node and released node algorithms and estimated the Wigner 
transition [26] at r s ss 100 ± 20 and a polarized liquid state stable for 75 < r s < 100. 
Calculations of Ortiz [22] found a lower value for Wigner crystallization of r s = 65 ± 10 
close to the older VMC results, in addition to a higher density spin-polarized phase 
starting at r s = 25 ± 5. Recent improved backflow calculations do not support these 
values. More accurate crystal calculations [32] estimate the melting transition at r s — 
106 ± 1 and the spin polarization transition [30] at r s = 50 ± 2. Melting of the bosonic 
Wigner crystal at non-zero temperature has been studied using Path Integral Monte 
Carlo [47] finding an abrupt transition between classical thermal melting and quantum 
mechanical pressure melting. 

Concerning the energies, a major improvement in the fluid phase was the introduction 
of backflow functions [19]. Twist averaging and analytic backflow[21] in the last few years 
has improved accuracy by potentially by an order of magnitude. We anticipate having 
much more accurate energies for the electron gas in the next few years. 

Correlation energy, structure factor, radial distribution function, and momentum dis- 
tribution of the spin-polarized uniform electron gas are tabulated in rcf. [33]. Analytic 
static structure factors and pair-correlation functions obtained by fitting the QMC re- 
sults and exact analytical limits for the unpolarized homogeneous electron gas including 
spin-resolved (up-up and up-down) correlation holes and static structure factors (unpo- 
larized, for the range 0.8 < rs < 10 are tabulated in ref. [34]. Quantum Monte Carlo 
calculations of the energy of the relativistic electron gas are considered in ref. [35]. A 
recent paper [36] considered the usual properties of the kinetic energy of the electron 
gas at finite temperatures. This followed observations reported in a restricted PIMC 



Introduction to Quantum Monte Carlo Methods Applied to the Electron Gas 



21 



calculations of the electron fluid [37]. 

For the 2DEG, the original VMC calculation found Wigner Crystallization[ll] at 
r s = 33 with a spin polarized state stable for 13 < r s < 33. This was followed up by 
a more accurate DMC study [38] resulting in an estimated transition at r s = 37 ± 5. 
Backflow trial functions were introduced [18]. There have been more recent DMC studies 
of the transition region and spin polarization [72, 74] with the conclusion that the spin 
polarized state may be stable between 26 < r s 35 though the energy differences are smaller 
than the systematic errors. We discuss next the calculation of Fermi liquid parameters 
in the 2DEG. Further studies of the correlation energy and spin polarization in the 2D 
electron gas are reported in ref. [39]. 

4'1. Fermi liquid parameters. - Now I briefly sketch the calculation of the energy 
of particle-hole excitations. We calculated the ground state and the lowest particle- 
hole excitations of the systcm[20] using correlated sampling techniques similar to that 
used in optimizing the trial wave function. The ground state consists of filled shells 
of the wave vectors allowed by periodic boundary conditions. To calculate the fermi 
liquid parameters, we consider excited states which are generated by exciting a single 
electron from the last occupied shell of the ground state to the first unoccupied shell 
(see Fig. 1 of Ref. [20] ) . Because these states have different total momentum, they will 
be orthogonal. Two different excitations, spin-parallel excitations and spin-antiparallel 
ones, are possible. It is efficient to use a guiding or sampling function of the form: 

(39) **,(R) = ao ^2 (R) + ^|^ (R) |2 5 

a 

where is the trial function for the ground state, ^ a for excited state a. The constant 
a is set to be roughly equal to the number of excitations considered. This guiding 
function is non-negative and zero only where all states under consideration have zeroes. 
One then uses a single very long random walk to calculate the energies E Q ,{E a } using 
the reweighting method discussed earlier. Then the energy differences will have a much 
smaller error, than do the individual energies. See ref. [20] for details. 

Following Fermi liquid theory analysis, the energy difference between two excited 
states a and (j is given by 

(40) AE a0 = Ui ± fi) [- «>s(Z0 a ) + cos(Z^)] , 

i=i 

where 9 a (p) is the angle between particle momentum k p and hole momentum k; t in 
excitation a{(3) and the +(— ) sign corresponds to parallel (antiparallel) spins between 
particle and hole. The effective mass is determined by the first-order spin-symmetric 
component /*: 

(41) £ = a - r^/rr 1 • 
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The ff can be obtained by 

(42) Nft 



NAEjl + NAEjj 
2(— cos 9i + cos 6*4) 



where AE^ 11 ' is the energy difference between spin-parallel (-antiparallel) excitations 
1 and 4. See ref. [20]for results. 

4'2. Polarization of electron gas. - In this section we discuss recent calculations of the 
polarization of the electron gas [30]. Despite years of active research, the properties of 
thermodynamic phases of the electron gas are not resolved experimentally or theoretically 
at intermediate densities. There has been recent interest in the low density phases 
spurred by the observation of a ferromagnetic state in calcium hexaboride (CaB 6 ) doped 
with lanthium[40]. The magnetic moment corresponds to roughly 10% of the doping 
density. The temperatures (600K) and densities (7 x 10 19 /cm 3 ) of this transition are 
in rough agreement with the predicted transition in the homogeneous electron gas. [22]. 
However, to make a detailed comparison, it is necessary to correct for band effects. For 
example, conduction electrons are located at the X-point of the cubic band structure and 
thus have a six-fold degeneracy. The effective mass of electrons at this point and the 
dielectric constant are also changed significantly from their vacuum values. These effects 
cast doubt on the viability of the electron gas model to explain the observed phenomena. 
Instead, excitonic and other models have been proposed to explain the ferromagnetism. 

Considering now the spin degrees of freedom, at small r s , electrons fill the Fermi 
sea with equal number of up spin and down spin electrons to minimize the total kinetic 
energy and thus the total energy; the system is in the paramagnetic state. As the density 
decreases and before the freezing transition, there is a possibility that the electrons 
become partially or totally polarized (ferromagnetic). The spin polarization is defined 
as £ = \Nf — ATj_ | /N, where and N± are the number of up and down spin electrons, 
respectively and N = Nf + N±. For paramagnetic phase £ = and for ferromagnetic 
phase ( = 1. 

This polarization transition was suggested by Bloch[41] who studied the polarized 
electronic state within the Hartree-Fock (HF) approximation. He found the ferromagnetic 
state favored over paramagnetic state for r s > 5.45, almost within the density of electrons 
in metals. However, HF is not reliable for r s > 0. 

Ceperley[ll] using VMC with a Slater-Jastrow trial function determined that the 
transition between the polarized and unpolarized phase occurred at r s = 26 ± 5. Using 
DMC[26] it was estimated that the polarized fluid phase is stable at r s — 75 ± 5. An 
extension to this work[42] found the £ = 0.5 partially polarized fluid becomes stable at 
roughly r s « 20 and the completely polarized state never stable. More recently Ortiz 
et al.[22] applied similar methods[22] to much larger systems (N < 1930) in order to 
reduce the finite-size error. They concluded the transition from the paramagnetic to 
ferromagnetic transition is a continuous transition, occurring over the density range of 
20 ± 5 < r s < 40 ± 5, with a fully polarized state at r s > 40. 
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Fig. 4. - Energy versus spin polarization at r a — 50 for 54 electrons using TA with 10 3 twist val- 
ues. Compared are calculations with SJ and BF-3B wavefunctions and with two QMC methods: 
VMC and DMC. 



Due to the very small energy differences between states with different polarizations, 
systematic errors greatly affect the QMC results, accurate trial wavefunctions are crucial 
to compute the small energy differences between different polarization states. Figure 4'2 
shows the energy vs. polarization at r s — 50 using different trial functions and simulation 
methods. The SJ trial function with VMC has the highest energy for all polarizations 
and at this level of accuracy finds the fully polarized phase to be stable, in agreement 
with earlier VMC calculations. [11] However, using the best BF-3B trial function, the 
variational energies are lowered significantly with the unpolarizcd energy dropping more 
than the polarized case so that the polarized phase is no longer stable. DMC calculations 
confirm this result. Note that the DMC energies determined using the NI phases (or 
nodes) give energies lower than the BF-3B variational energies confirming the importance 
of accurate DMC calculations. The use of BF-3B wavefunctions with DMC leads to the 
lowest ground state energies, hopefully, very close to the exact energy. Twist averaging 
is particularly advantageous for polarization calculations since shell effects dominate the 
polarization energy. 
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Fig. 5. - The spin polarization energy of the 3DEG times r 3 s ^ 2 in Ry/electron at various densities 
using a cubic polynomial fit to the QMC data. The density, r s , is denoted on the right axis. 



The results are shown in figure 5. A polarization transition is evident. At r s = 40, the 
system is still paramagnetic, with the unpolarized phase stable. As the density decreases, 
at r s w 50, the system becomes unstable with respect to spin fluctuations. The partially 
polarized states become stable at r s > 60. As the electronic density continues to decrease, 
the fully polarized state has a lower energy with respect to unpolarized state at r s > 80, 
however, we find that the partially polarized state has an even lower energy. 

In Fig. 6 is shown the predicted square of the optimal polarization versus density. 
It is found that the equilibrium polarization is described by £ 2 = (r s — r* ) /62 with the 
critical density r* = 50 ± 2. As the density decreases, the stable state becomes more 
and more polarized, approaching fully polarized at freezing, r s ss 100. Quantum critical 
fluctuations, not present in systems with N < 162, could modify the behavior of the spin 
polarization energy near the critical density. 

The quoted error bar on the critical density estimates the statistical errors, not the 
systematic errors arising from the fixed-phase approximation. The experimental and 
theoretical results on polarized helium give caution on placing too much confidence in 
the estimate of the polarization transition. Even using the accurate optimized BF-3B 
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Fig. 6. - The square of the spin polarization versus r s . The curves were obtained using fits 
in the Fig. (4). The line is a fit though the points. The value at r s = 40 was obtained by 
extrapolation from physical values of £. 



wavefunctions, the magnetic susceptibility in liquid 3 He does not agree with experiment 
at low pressure and the polarized phase is nearly degenerate with the unpolarized phase 
at the freezing density. [75] The present results also do not preclude the existence of phases 
with other order parameters such as superffuids, as occurs in ground state of liquid 3 He. 
In fact, it is rather likely that the ground state of the electron gas will have such a phase 
at the lowest fluid densities. 

One can use the calculated energies to estimate the finite temperature behavior within 
the Stoner model[44] arising in the theory of itinerant magnetism [45]. The Stoner model 
differs from HF by replacing the Coulomb interaction by a zero range one, a repulsive 
delta function potential: £\ . gS(rij). One can view this approach as the first step to a 
full Fermi-liquid description of the quasiparticle interactions, and use the QMC data to 
determine the strength of those interactions. One expects that the bosonic correlations 
have screened off the long-range interaction, leaving only a short-ranged spin-dependent 
term that can be modelled by a contact interaction. 

In the Stoner model, the energy is evaluated within the mean field approximation. 



26 



D. M. Ceperley 



The energy at zero temperature in the thermodynamic limit is: 

(43) E oc (1 + C) 5/3 + (1 - C) 5/3 + 0.054gr s 2 (l - ( 2 ). 

For gr 2 s < 20.5 the system has an unpolarized ground state and for gr 2 > 24.4 the ground 
state is ferromagnetic. For intermediate couplings, the ground state has a partial spin 
polarization at zero temperature, similar to the observed behavior of the electron gas at 
low density. 

Although the polarizations are qualitatively correct, the above functional form does 
not fit well the DMC data. In addition, assuming g docs not have a very strong density 
dependence, the Stoner model predicts that the partially polarized density range should 
be quite narrow, from 50 < r s < 54, while as the QMC results indicate a much broader 
density range. Certainly, the assumption of a zero-range interaction of quasiparticles 
is too restrictive. However, we note that in the case of the 3D Ising model, the mean 
field estimate of the critical temperature is approximately 20% greater than the exact 
value, suggesting that the Stoner model will give a reasonable estimate of the transi- 
tion temperature if the effective couplings are determined from the QMC ground state 
energies. 

We use the Stoner model to make a estimate of the transition temperature of the 
polarized phase as follows. The free energy [46] in a fixed volume V in the Stoner model 
is: 

(44) F = Fo (N l ) + F (N r ) + ^^. 
For large N the free particle free energy is: 

(45) F (N) = Nfi - k B T^2ln(l + e- fi(ek -^). 

k 

The chemical potential /i of each spin species is determined by the number of particles 
with that spin. 

Fig. 7 is the estimated phase diagram of the 3D electron gas. Note that both the tem- 
peratures and densities of our calculated magnetic transition are four orders of magnitude 
smaller than that found experimentally [40] in CclBq. Even assuming errors because of 
uncertainties in material properties and systematic errors in the calculation of T c , these 
estimates are very difficult to reconcile with experiment, apparently ruling out an elec- 
tron gas model of fcrromagnctism in this material. An estimate of the limit of stability 
of the Wigner crystal [47] is also shown. 

There is recent work on the low density electron gas in 2D, studying both the Fermi 
liquid [39] and the Wigner crystal phase. Using QMC similar techniques in the liquid 
phase, a polarized phase is found to be stable between 26 < r s < 35, though the energy 
differences are even smaller than in 3D. The partially polarized phase is never stable. 
In the 2D Wigner crystal, path integral methods[49] were used to derive directly the 
spin Hamiltonian. It was found that the ground magnetic state is a spin liquid though 
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Fig. 7. - The predicted phase diagram of the 3DEG from QMC calculations. The solid line 
the mean-field estimate of the magnetic transition temperature from the Stoner model, where 
the spin interaction is estimated from the zero temperature QMC data. The dotted line is the 
energy difference between the unpolarized and partially polarized system. 



the ferromagnetic state has only a slightly higher energy at melting. We will discuss 
this below. Analogous calculations of the magnetic phase diagram of the WC in 3D are 
underway [50]. 

4'3. Calculation of spin and charge response. - The majority of QMC calculations have 
been for equilibrium properties such as energy, one-particle-orbital occupation numbers, 
and static correlation functions. Here I briefly discuss the basis of calculations of the 
spin and charged response[52, 53]. Apart from their intrinsic interest, these response 
functions are of importance to density functional developments beyond LDA[54, 55]. 

The static density-density response function is directly calculable by QMC. Rather 
than evaluate it in terms of the fluctuation-dissipation theorem, we use the definition of 
static response function by applying a static external potential: 



(46) 



v e xt(r) = 2vqCos(q ■ r), 
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which induces a modulation of the density with respect to its mean value, n . Such 
a modulation contains periodic components at all wavevectors that are non-vanishing 
integer multiples of q. In particular, one finds a modulation with wavevector q, ni(r) = 
2n q cos(q • r), where 

(47) n q = x(q)vq + C 3 vl + • • • , 

only contains odd powers of u q . Here x{l) denotes the static density-density linear 
response function in Fourier space. Similarly the ground state energy (per particle) can 
be expanded in even powers of w q : 

(48) E v = E +^vl + C 4 v* + ---. 

The coefficients C 3 and C4 in the above equations are determined by the cubic response 
function. QMC allows the direct evaluation of both n q and E v , for given q and w q . One 
performs simulations at few coupling strengths w q and then extract x(<l) as we U as the 
higher order response functions from the calculated n q or E v , by fitting in powers of w q . 

There are several possibilities for the trial function in the external field. The first 
form is: 



(49) *r(B-) = *t(R) II ex P\ucos{ci • r,)], 

i 

with a a new variational parameter, related to the external potential strength t> q . One 
can easily show that ^y(R), to leading order in a, correctly yields n(r) = no + 2ajcos(q- 
r), with 7 a function of the density uq. To determine the relationship between the 
amplitude of the external potential u q , and the variational wavefunction parameter a, 
it is convenient to fix the wavefunction parameter a and determine the corresponding 
potential w q which satisfies the minimum condition dE v /da = 0. Once the trial function 
has been optimized, DMC is used to evaluate the total energy and the Fourier 
component n q of the density. One should keep in mind that while the estimate of 
the energy is exact, within statistical errors, the extrapolated estimate [15] yielding n q is 
approximate. In fact its accuracy depends quadratically on the difference 5^ = ^0 — ^t- 
Therefore, it is better to evaluate x(g) from the energy, using Eq.(48). 

Although the nodal structure of is thought to be accurate enough for the un- 
perturbed system, the nodes of the trial function discussed above do not depend on the 
perturbation; one knows is incorrect for non-interacting particles. A second trial function 

(50) <&£ ;2 = D]D\e- u , 

is obtained by constructing the determinants D v s in terms of one-particle orbitals (Math- 
icu functions) for non-interacting electrons in a new external field, v'(r) = acos(q ■ r), 
with a a variational parameter. In general, <J>y 2 wm possess a different nodal structure 
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, while still yielding a modulated density. It is found that the energy of this second trial 
function is lower. Note that the effective potential should be weak enough that the filling 
of the single particle orbitals is not altered. 

This method of evaluating the static response functions is computationally demand- 
ing: to obtain x(g), at a given thermodynamic state, requires a few simulations for each 
wavevector q — thus involving of the order of tens of simulations to construct over 
the range of relevant wavevectors. Though the use of the fluctuation-dissipation theorem 
could in principle require less computer time, since the entire response function would 
be calculated at once, estimation of the systematic errors for the convergence of the time 
integral arc problematical and the formula is not appropriate with fixed-node DMC. See 
refs. [52, 53] for details and results of calculations of the response functions. 

5. Path Integral Monte Carlo and Applications 

Up to this point we have considered zero temperature methods, variational Monte 
Carlo and Diffusion Monte Carlo. Path Integrals have some significant advantages over 
zero temperature methods, as well as disadvantages, of course. Among the advantages 
is the absence of a trial wavefunction which means that quantum expectation values, 
including ones not involving the energy, can be directly computed. For the expert, the 
lack of an importance function may seem a disadvantage; without it one cannot push 
the simulation in a preferred direction. However, as the quantum system becomes more 
complex, it becomes increasingly difficult to devise satisfactory trial functions. It is my 
opinion that if quantum simulations are to attack the type of complex physical situations 
that classical simulations routinely deal with, it is better to have a formulation without 
a trial function. Only the Hamiltonian should enter. Of course, the explicit formulation 
at finite temperature also makes comparison with experiment more direct. 

Thermodynamic properties are averages over the thermal TV-body density matrix 
which is defined as a thermal occupation of the exact eigenstates <j>i(R): 



(51) 




The partition function is the trace of the density matrix. 



(52) 




Other thermodynamic averages are obtained as: 



(53) 




The density matrix can be calculated using path integrals. As first shown by Fcynman 
(1953), the many-body density matrix can be obtained using classical-statistical methods 
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on polymer-like systems. The density matrix of a many-body system at a temperature 
k B T — can be written as an integral over all paths {Rt}- 



The path R(t) begins at VRq and ends at Rp, and V is a permutation of particle labels. 
For N particles, the path is in 3iV dimensional space: R t — (i"i t ,r 2 t . . . rjvt)- The upper 
sign is to be used for bosons and the lower sign for fermions. The action of the path, 
S[R t ], is given by: 



Thermodynamic properties, such as the energy, are related to the diagonal part of the 
density matrix, so that the path returns to its starting place or to a permutation V of 
its starting place after a "time" (3. 

Since the imaginary-time action S[R t ] is a real function of the path, for boltzmannons 
and bosons the integrand is nonnegative. It can thus be interpreted as a probability 
of an equivalent classical system and the action as the classical potential energy of a 
"polymer." To perform Monte Carlo calculations of the integrand, one makes imaginary 
time discrete, so that one has a finite (and hopefully small) number of time slices and thus 
a classical system of N particles in M time slices; an equivalent NM particle classical 
system of "polymers." If the path integral is performed by a simulation method, such as 
a generalization of Metropolis Monte Carlo or with Molecular Dynamics, one can obtain 
essentially exact results for systems such as the properties of liquid 4 He at temperatures 
near the superfluid phase transition, the exchange frequency in quantum crystals, and 
quantum particles immersed in classical systems [61]. 

Note that in addition to sampling the path, the permutation is also sampled. This is 
equivalent to allowing the ring polymers to connect in different ways. This macroscopic 
"percolation" of the polymers is directly related to superfluidity as Feynman (1953) first 
showed. Any permutation can be broken into permutation cycles, i. e. into 2-, 3-, ... 
exchange cycles. Superfluid behavior can occur at low temperature when the probability 
of exchange cycles on the order of the system size is nonncgligiblc. For more details on 
the path integral theory of Bose superfluids and how one carries out the Monte Carlo 
calculations see[61]. 

However, the straightforward application of those techniques to Fermi systems means 
that odd permutations subtract from the integrand. This is the "fermion sign problem." 
Path integral methods as rigorous and successful as those for boson systems are not yet 
known for fermion systems in spite of the activities of many scientists throughout the 
last four decades. 

Now let us consider how particle statistics are expressed in path integrals. For sys- 
tems of identical particles, the states can be classified into symmetric and antisymmetric 



(54) 




(55) 
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states. The fermion density matrix is denned by restricting the sum to be only over 
antisymmetric states. (Similarly for other symmetries such as momentum or spin.) We 
shall denote the statistics of the particles by subscripts: p F will denote the fermion den- 
sity matrix, p B the boson density matrix, pu the boltzmannon (distinguishable particle) 
density matrix, and p any of the above density matrices. 

For the moment, let us consider a single component system of spinlcss fermions. (By 
spinlcss fermions, we mean that the spatial wavefunction is antisymmetric with respect 
to interchange of all pairs of spatial variables.) Let V be one of the Nl permutations of 
particle labels. Using Eqs. (51) the density matrix has the following symmetries: 

p(R,R';/3)=p(R',R;/3) 

(56) p F (R, R'; (5) = (-lfp F (VR, R'; 0) 

p F (R,R';(3) = (-lfp F (R,VR';f3). 

One can use the permutation (or relabeling) operator to construct the path integral 
expression for the boson or fermion density matrix in terms of the Boltzmann density 
matrix: 

(57) p B/F (R, R'; /?) = ^ £(±1) 7 #\ P)- 

' V 

More generally, one uses some projection operator to select a desired set of states from 
the distinguishable particle density matrix which contains all states. From Eq. (56) we 
could (anti)symmetrize with respect to the first argument, the last argument or both. 
This connection between the boltzmannon density matrix and the Bose/Fermion density 
matrix is important because it is the boltzmannon density matrix that is built naturally 
from paths. 

An alternative definition of the density matrix is by its evolution in imaginary time, 
the Bloch equation: 

(58) _ d P {R,R';t) =Hp{RR ,. t) 

which obeys the boundary condition at t = for boltzmannon statistics: 

(59) p D (R,R';0) = 5{R- R') 
or for Bose or Fermi statistics: 

(60) PB/F (R, R'; 0) = ^ ^(±lf6(PR - R'). 

' V 

The high temperature boundary condition is an (anti)symmetrized delta function. 
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Path integrals are constructed using the product property of density matrices: 

(61) p{R ,R 2 ,Pi + (3 2 ) = j dR lP {R ,R 1 ;(3 l )p{R 1 ,R 2 ;l3 2 ). 

The product property holds for any sort of density matrix. 

If the product property is used M times we can relate the density matrix at a tem- 
perature to the density matrix at Mf}^ 1 : 

(62) Pb/f(Ro, Rm, 0) = ]r } H(± 1 ) 7 ' fdRt... dR M -i 

■ v 

Pd(VR ,Ri;t) . ..pd(R M -i,Rm;t). 
The sequence of intermediate points {Ri, R 2 , ■ ■ ■ , Rm-i} is the path, and the time step 

is T = 0/M. 

As the time step gets sufficiently small, we can write down an explicit expression 
for the density matrix po and thereby an explicit expression for p(R , Rm', ft), but one 
with lots of intermediate integrals and the permutational sum to perform. The Trotter 
theorem tells us that for sufficiently small r we can assume that the kinetic and potential 
operators commute so that: e~ T ^ = e Tl e~ rV '. Define the (boltzmannon) action as 
Sd(R, R'; t) = — ln[p D (R, R'; r)]. Then for small time step the action is: 

(63) S D (R, RI; r) = ^ ln^Ar) + {R ^ ? + \{V(R) + V(R')). 

This is known as the primitive approximation to the action. The form of the action is 
analogous to the potential energy of classical "polymer" system with harmonic springs 
between neighboring beads and an interpolymer potential between different chains. 

The boson action is real, but it is expressed as a sum over permutations. For large 
N it is not possible to evaluate directly the sum since it has Nl terms. It is better to 
leave the bosonic symmetrization as an explicit boundary condition on the paths and to 
sample the permutations as well as the paths. We will follow the same philosophy with 
restricted fermion paths, the reasons not being the difficulty of evaluating the resulting 
determinant (that is easier for fermions than for bosons) but to avoid the minus signs. 
For bosons or fermions we can (anti)symmctrize anywhere along the path as many times 
as we like. However relabeling is only necessary once. By convention we will just relabel 
the first step. 

In the direct fermion method one sums over permutations just as for bosonic systems. 
Odd permutations then contribute with a negative weight. The direct method has a 
major problem because of the cancellation of positive and negative permutations as noted 
by Feynman and Hibbs (1965). It is possible to demonstrate that a distribution having 
both positive and negative regions will have an exponentially vanishing signal/noise ratio 
in a Monte Carlo calculation. At high temperatures 



(64) 



£ = cxp[-2p7V(27rA/3) 3/2 ]. 
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At the degeneracy temperature: £ w exp(— N). The direct fermion method, while exact, 
becomes exceedingly inefficient as (3 and N increase, precisely when the physics becomes 
interesting. 

The restricted path identity (or fixed-node) shows that the nodes of the exact density 
matrix determine the rule by which one can take only paths with the same sign; we can 
arrange things so that we only get positive contributions. We only use paths which do not 
cross the nodes of the fermion density matrix, and the results are identical to allowing 
all of the paths. Though the identity is exact, one faces precisely the same difficulty 
as in zero temperature methods: the nodal surfaces of the interacting fermion density 
matrix are unknown, so nodal surfaces coming from cither perturbation theory, or the 
variational principle must be used. These methods lead to a very interesting description 
of fermion systems in terms of exchanging paths. However, we do not have space to 
describe them here and refer instead to the litcrature[57, 58] 

5T. Defects in the 2D Wigner Crystal. - At sufficiently low densities, the ground 
state of a 2D electron gas is a Wigner crystal. In this section, I summarize the results 
of calculations of point-like defects in a clean 2D electron Wigner crystal [62]. The moti- 
vation for the study of defects is two- fold. First, localized defects are present in a finite 
concentration at any nonzero temperature (they have even been speculated to exist at 
zero temperature a super solid). Secondly, the melting process in 2D can be influenced 
or even determined by defect formation. 

At large values of r s , the exchange contributions to the energy are small as we discuss 
below, so one can assume distinguishable electrons and neglect antisymmetry but include 
it for r s < 75. However, to keep the system stable, one must forbid particle exchange 
by enforcing the condition |rj — Sj| < 1.1a where a is the nearest neighbor distance 
and Si is the i th lattice site. Such a "tether" is realistic because in a quantum crystal 
exchanges are rare and the wave function is peaked around the lattice sites. One can 
also use restricted path integrals[57] to account for Fermi statistics. Only paths entirely 
in the positive region of the Slater determinant are allowed. Such a restriction is exact 
if the nodal surfaces of the trial density matrix are correct. For the nodes one can use 
a Slater determinant of Gaussian orbitals exp(— c(r^ — Sj) 2 ) with Sj the lattice sites and 
c taken from [11] with a ferromagnetic spin arrangement. Calculation of the tunnelling 
frequencies to determine the ground state magnetic ordering indicate that the magnetic 
energies are always much less than the defect energies [49] and that the system is nearly 
ferromagnetic in near melting. We find that the restriction also serves to stabilize the 
crystal against melting. 

The energy to create N de { defects in a system with N t = N lattice sites is 

(65) AE dcf = [e(N + N dcf ) - e(N)] (N + N dei ), 

where e(n) is the energy per electron for a system containing n electrons with area 
A = n/p. This was evaluated with independent PIMC calculations with different number 
of particles rather than more efficient procedures where particles are inserted or removed. 
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Fig. 8. - a) Hexagonal 2D Wigner crystal with a single 6-coordinated vacancy and a 3-fold 
coordinated centered interstitial (CI) defect, b) Pair vacancy and CI defects. 



Such differential procedures are difficult because of the combination of large relaxation of 
the lattice, the very large zero point motion and the antisymmetry. The calculations are 
time-consuming as the system gets larger because one needs high accuracy to obtain the 
energy difference but the direct method allows better control over the systematic error. 

Fig. 8 shows a vacancy and centered-interstitial (CI) defects in a 2D hexagonal lattice. 
Fig. 9 demonstrates that the formation energy for an interstitial is consistently lower 
than creation energy for vacancies in 2D Wigner crystals at all densities. For a Coulomb 
system without fermion exchange, the energy of a defect can be expanded as 

(66) E D = ar- 1 + c m r-V 2 + c 2 r- 2 . . . 

The first term is the static potential energy of the defect, the second the harmonic 
energy of the defect. Cockayne and Elser[59] have done exact calculation of ci and c 3 / 2 . 
Shown in the lower panel of Fig. 9 is the anharmonic contribution defined as the excess 
energy beyond the harmonic calculation. Anharmonic effects lower the defect energies 
to approximately half the harmonic value at r s = 50. Cockayne and Elser find that the 
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Fig. 9. - Formation energy (top figure) for vacancy (circle) and CI (square) defects as a function 
of l/r s at a temperature of 1.25xlO _5 Ry for a system of 120 lattice sites (open and solid 
syimbols are for Boltzman and Fermi statistics respectively). The triangles are results with 
340 lattice sites. The lines (solid for vacancy and dashed for interstitial) are from harmonic 
calculations [59]. The anharmonic energy is shown on the bottom panel. 



defect energy vanishes for CI defects at r s = 15 ± 1 and r s = 9 ± 1 for vacancies. We 
estimate that Ejj vanishes for intcrstitials at r s = 35 ± 2 and r s = 29 ± 2 for vacancies. 
The proximity of the vanishing of the interstitial creation energy to the melting density 
is highly suggestive that interstitial defects play a role in the melting process. The 
proliferation of vacancies for r s < 40 could result in a continuous rather than first order 
quantum melting, analogous to what happens for the classical 2D Wigner crystal. It is 
possible that the assumption of the ferromagnetic spin arrangement has stabilized the 
crystal with respect to interstitials for 35 < r s < 40. 

5'2. PIMC of exchange frequencies and magnetic ordering in the 2 wigner crystal. - 
In this section I briefly review calculations [49] to determine the spin-spin interaction in 
the Wigner crystal, using Thouless'[63] theory of exchange. According to this theory, 
in the absence of point defects, at low temperatures the spins will be governed by a 
Hamiltonian of the form: 

(67) H S pin = - y~^(— 1) P JpVspin 

p 



36 



D. M. Ceperley 



2d 



10- 3 








-s 


10- 4 
10- 5 

10" 6 


I / 
± 




quantum crystal--''" ^ 


quanturfi 
fluid _; 


10- 7 










10- 8 






,-* / 

/ / 
- - / 




10- 9 






// spin liq 




10" 10 




/ 






10- 11 




/. 

/ \ 
il ' 






io- 12 




II 






10- 13 




il 

h 






1 Q-14 




;h 1 







12 3 

100/r 



Fig. 10. - Phase diagram of the 2dEG. The estimated melting line is based on Lindemann's 
criteria. The long-dashed line represents the Debye temperature. The dotted line is the Curie- 
Weiss constant 0, the short dashed line is coefficient of the specific heat at high temperature, 
J c as defined in the text. The vertical line is the estimated zero temperature ferromagnetic (F) 
transition. 



where the sum is over all cyclic (ring) exchanges described by a cyclic permutation P, Jp 
is its exchange frequency and V sp i n is the corresponding spin exchange operator. Path 
Integral Monte Carlo (PIMC) as suggested by Thouless[63] and Roger[64] has proved 
to be the only reliable way to calculate these parameters. The theory and computa- 
tional method have been tested thoroughly on the magnetic properties of bulk helium 
obtaining excellent agreement with measured properties [61]. Rather surprisingly, it has 
been found [65] that in both 2D and 3D solid 3 He, exchanges of 2, 3 and 4 particles have 
roughly the same order of magnitude and must all be taken into account. This is known 
as the multiple spin exchange model(MSE). 

A WKB calculation of the exchange frequencies in the 2dWC by Roger [64] predicted 
that the three electron J3 nearest neighbor exchange would dominate, leading to a fer- 
romagnetic^) ground state. Recent calculations [66, 67] have confirmed and extended 
those of Roger. 

An exact method for calculating the exchange frequency in quantum crystals has been 
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Fig. 11. - Exchange frequencies versus r x J 2 . The solid line is 10 3 of the kinetic energy. 



previously developed and applied to solid 3 He [68, 61, 69]. One computes the free energy 
necessary to make an exchange beginning with one arrangement of particles to lattice 
sites Z and ending on a permuted arrangement PZ: 

(68) F P (/3) = Q(P,l3)/Q(I,l3)=Unh(Ml3-(3 )). 

Here Q(P,f3) is the partition function corresponding to an exchange P at a temperature 
1//3. I is the identity permutation. Note that these paths are of "distinguishable" 
particles since Fermi statistics are implemented through the spin Hamilitonian in Eq. 
(1). The function Fp(f3) is found with the Bennett method, which directly calculates 
free energy differences and thereby determines Jp and /3o- 

The exchange frequencies vary rapidly with density as shown in fig. 11. One can see 
that they are much less than the zero point energy of the electrons, thus justifying the use 
of Thouless' theory. The WKB method [65], where one approximates the Path Integral 
in Eq. (2) by the single most probable path, explains most of this density dependence. 
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In the 2dWC, the WKB expression for the exchange frequency[64, 66, 67] is: 
(69) J P = Apir^r^e-^ 2 . 

1/2 

Here bpr s is the minimum value of the action integral along the path connecting PZ 
with Z. The 3 particle exchange exponent is the smallest indicating that as r s — > oo, 
J 3 will dominate and the system will have a ferromagnetic ground state. However, note 
that in Fig. 11 J 2 > J 3 for r s < 90. 

Having determined the exchange frequency, one is left with the spin Hamiltonian of 
Eq. (1). This is a non-trivial many-body problem. For spin 1/2 systems, J 2 an d J3 
contribute only with a nearest neighbor Heisenberg term: J| fl = J2 — 2J 3 . This term 
is negative (ferromagnetic) but approaches zero near melting. For convenience we use 
J4 as a reference to fix the overall scale of the magnetic energy Neglecting J &p , the 
Hamiltonian has 3 remaining parameters J% s /J4, J5/J4 and Jg/i/ Ji- The dependence of 
these ratios on density is shown in fig. (5). 

High temperature series expansions [70] determine the the specific heat CV and mag- 
netic susceptibility Xo/x f° r temperatures ksT 3> Jp- The susceptibility is given by: 
Xo/X = T -6 + B/T ... and the specific heat C v /Nk B = (3J C /2T) 2 + . . . where the 
Curie-Weiss constant is given by 9 = — 3(J| ff + 3J4 — 5J5 + b/SJ^h + 15/4J 6p ) with a 
quadratic expression of the J's for J c . These two constants, which set the scale of the 
temperature where exchange is important, are shown as a dotted and dashed lines on fig. 
(1). Note that 9 changes from positive to negative at r s w 130. Both 9 and J c decrease 
very rapidly at low density showing that experiments must be done at r s < 60 if spin 
effects are to be at a reasonable temperature, e.g. T c > O.lmK (assuming the values for 
Si-MOSFET: e = 7.7, m* = 0.2). 

The zero temperature state can be studied by exact diagonalization[71] (ED) of an 
N site system. (The present limitation is N < 36.) Two phases are important: the 
ferromagnetic (F) phase and an antiferromagnetic(AF) phase. The F-AF transition is 
shown in fig. 10. The ferromagnetic phase is obtained only at very low density: the 
transition for the 2dWC will occur at r s = 175 ± 10. At higher density, the frustration 
between large cyclic exchanges (4-6 body loops) results in a disordered spin state [71]. 
For example, at r s — 100, is a spin liquid with a gap to all excitations. At higher 
densities (r s < 100), the trajectory of the MSE models parallels the F-AF phase line, 
with the possibility of a re-entrant ferromagnetic phase for r s < 40. Note that QMC 
calculations [72] of the normal fcrmi liquid at r s = 30 show that the ferromagnetic phase 
has a slightly lower energy that the unpolarized phase. Hence both the high density 
2dWC and the low density electron fluid are characterized by a spin Hamiltonian which 
is nearly ferromagnetic. 

if if if 
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with the facilities at the NCSA. This article summarizes a series of calculations and 
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Fig. 12. - Spin Phase diagram of the 2DEG as a function of exchange ratios. The dotted line is 
the flow of spin hamiltonian space versus r a (lower numbers); also shown the are estimated values 
of J&h/ Ji (upper numbers). The solid lines are the limit of the ferromagnetic phase according to 
ED [71] at J 6h /J4. = {0.2,0.4,0.6}. The 2dWC crosses into the F region for r s w 175. The (*) 
are empirical estimates of the spin Hamiltonian of 2d 3 He at several densities [73]. 
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